I am largely following the DADA2 tutorial by Benjamin Callahan that can be found here.
## [1] "A1" "A10" "A11" "A12" "A2" "A3" "A4" "A5" "A6" "A7" "A8" "A9"
## [13] "B1" "B10" "B11" "B12" "B2" "B3" "B4" "B5" "B6" "B7" "B8" "B9"
## [25] "C1" "C10" "C11" "C12" "C2" "C3" "C4" "C5" "C6" "C7" "C8" "C9"
## [37] "G5" "G6"
Viewing the quality of 2 of the forward samples
Viewing the quality of the same 2 of the reverse samples
You can see in above figures that the quality of reads drop off towards the end, so we need to filter out these low quality reads
## reads.in reads.out
## A1_16S_R1.fastq 46968 44461
## A10_16S_R1.fastq 25343 23991
## A11_16S_R1.fastq 27398 26006
## A12_16S_R1.fastq 45822 44152
## A2_16S_R1.fastq 99390 94836
## A3_16S_R1.fastq 77568 73174
After filtering out the low quality reads, we have maintained about 95.5 of the original reads.
Viewing the quality of the same 2 of the forward samples post filtering and trimming
Viewing the quality of the same 2 of the reverse samples post filtering and trimming
The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. We want the estimated error rates (black line) to be a good fit to the observed rates (points), and the error rates to drop with increased quality.
We are now ready to apply the core sample inference algorithm to the filtered and trimmed sequence data.
Forward sample inference
## Sample 1 - 44461 reads in 9515 unique sequences.
## Sample 2 - 23991 reads in 6094 unique sequences.
## Sample 3 - 26006 reads in 6312 unique sequences.
## Sample 4 - 44152 reads in 7218 unique sequences.
## Sample 5 - 94836 reads in 13509 unique sequences.
## Sample 6 - 73174 reads in 9644 unique sequences.
## Sample 7 - 22069 reads in 4917 unique sequences.
## Sample 8 - 76841 reads in 25268 unique sequences.
## Sample 9 - 41897 reads in 7509 unique sequences.
## Sample 10 - 19892 reads in 5018 unique sequences.
## Sample 11 - 48504 reads in 7342 unique sequences.
## Sample 12 - 32048 reads in 7630 unique sequences.
## Sample 13 - 90475 reads in 7748 unique sequences.
## Sample 14 - 24910 reads in 5725 unique sequences.
## Sample 15 - 34782 reads in 5737 unique sequences.
## Sample 16 - 50069 reads in 5167 unique sequences.
## Sample 17 - 76155 reads in 7895 unique sequences.
## Sample 18 - 34374 reads in 5211 unique sequences.
## Sample 19 - 90460 reads in 8319 unique sequences.
## Sample 20 - 23566 reads in 5815 unique sequences.
## Sample 21 - 53651 reads in 5991 unique sequences.
## Sample 22 - 110660 reads in 10438 unique sequences.
## Sample 23 - 14928 reads in 3601 unique sequences.
## Sample 24 - 30116 reads in 6468 unique sequences.
## Sample 25 - 27927 reads in 5388 unique sequences.
## Sample 26 - 73356 reads in 7647 unique sequences.
## Sample 27 - 23594 reads in 6498 unique sequences.
## Sample 28 - 20343 reads in 4504 unique sequences.
## Sample 29 - 25590 reads in 6768 unique sequences.
## Sample 30 - 10144 reads in 3129 unique sequences.
## Sample 31 - 42054 reads in 7901 unique sequences.
## Sample 32 - 25161 reads in 7367 unique sequences.
## Sample 33 - 73293 reads in 8481 unique sequences.
## Sample 34 - 37175 reads in 5029 unique sequences.
## Sample 35 - 30216 reads in 4299 unique sequences.
## Sample 36 - 43608 reads in 5216 unique sequences.
## Sample 37 - 1795 reads in 421 unique sequences.
## Sample 38 - 162 reads in 63 unique sequences.
Inspecting the returned dada-class object for the first forward sample:
## dada-class: object describing DADA2 denoising results
## 144 sequence variants were inferred from 9515 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16
The DADA2 algorithm inferred 144 true sequence variants from the 9515 unique sequences in the first sample. There is much more to the dada-class return object than this (see help("dada-class") for some info), including multiple diagnostics about the quality of each denoised sequence variant, but that is beyond the scope of an introductory tutorial.
Reverse sample inference
## Sample 1 - 44461 reads in 9175 unique sequences.
## Sample 2 - 23991 reads in 5924 unique sequences.
## Sample 3 - 26006 reads in 6092 unique sequences.
## Sample 4 - 44152 reads in 7108 unique sequences.
## Sample 5 - 94836 reads in 12281 unique sequences.
## Sample 6 - 73174 reads in 8686 unique sequences.
## Sample 7 - 22069 reads in 4797 unique sequences.
## Sample 8 - 76841 reads in 24793 unique sequences.
## Sample 9 - 41897 reads in 7149 unique sequences.
## Sample 10 - 19892 reads in 4939 unique sequences.
## Sample 11 - 48504 reads in 7069 unique sequences.
## Sample 12 - 32048 reads in 7283 unique sequences.
## Sample 13 - 90475 reads in 7578 unique sequences.
## Sample 14 - 24910 reads in 5526 unique sequences.
## Sample 15 - 34782 reads in 5706 unique sequences.
## Sample 16 - 50069 reads in 5134 unique sequences.
## Sample 17 - 76155 reads in 7407 unique sequences.
## Sample 18 - 34374 reads in 4744 unique sequences.
## Sample 19 - 90460 reads in 8211 unique sequences.
## Sample 20 - 23566 reads in 5706 unique sequences.
## Sample 21 - 53651 reads in 5681 unique sequences.
## Sample 22 - 110660 reads in 10194 unique sequences.
## Sample 23 - 14928 reads in 3754 unique sequences.
## Sample 24 - 30116 reads in 6446 unique sequences.
## Sample 25 - 27927 reads in 5233 unique sequences.
## Sample 26 - 73356 reads in 7528 unique sequences.
## Sample 27 - 23594 reads in 6469 unique sequences.
## Sample 28 - 20343 reads in 4401 unique sequences.
## Sample 29 - 25590 reads in 6315 unique sequences.
## Sample 30 - 10144 reads in 3041 unique sequences.
## Sample 31 - 42054 reads in 7772 unique sequences.
## Sample 32 - 25161 reads in 7026 unique sequences.
## Sample 33 - 73293 reads in 8758 unique sequences.
## Sample 34 - 37175 reads in 5235 unique sequences.
## Sample 35 - 30216 reads in 4296 unique sequences.
## Sample 36 - 43608 reads in 5159 unique sequences.
## Sample 37 - 1795 reads in 465 unique sequences.
## Sample 38 - 162 reads in 69 unique sequences.
We now merge the forward and reverse reads together to obtain the full denoised sequences. Merging is performed by aligning the denoised forward reads with the reverse-complement of the corresponding denoised reverse reads, and then constructing the merged “contig” sequences. By default, merged sequences are only output if the forward and reverse reads overlap by at least 12 bases, and are identical to each other in the overlap region (but these conditions can be changed via function arguments).
The mergers object is a list of data.frames from each sample. Each data.frame contains the merged sequence, its abundance, and the indices of the forward and reverse sequence variants that were merged. Paired reads that did not exactly overlap were removed by mergePairs, further reducing spurious output.
We can now construct an amplicon sequence variant table (ASV) table, a higher-resolution version of the OTU table produced by traditional methods.
##
## 221 222 223 227 228 239 246 249 250 251 252 253 254 255 256 257
## 1 2 4 4 1 1 1 1 1 10 105 3259 140 9 5 2
The sequence table is a matrix with rows corresponding to (and named by) the samples, and columns corresponding to (and named by) the sequence variants. This table contains 3546 ASVs.
After viewing the distribution of read lengths, it looks like we have some that fall outside the expected range (244 - 264) so we will go ahead and remove these non target length sequences.
##
## 246 249 250 251 252 253 254 255 256 257
## 1 1 1 10 105 3259 140 9 5 2
This updated table now contains 3533 ASVs.
The core dada method corrects substitution and indel errors, but chimeras remain. Fortunately, the accuracy of sequence variants after denoising makes identifying chimeric ASVs simpler than when dealing with fuzzy OTUs. Chimeric sequences are identified if they can be exactly reconstructed by combining a left-segment and a right-segment from two more abundant “parent” sequences.
A total of 40 bimeras were identified from the 3493 input sequences, thus retaining 99.8% of sequences.
As a final check of our progress, we can look at the number of reads that made it through each step in the pipeline:
| Genotype | Treatment | Input | Filtered | Denoised Forward | Denoised Reverse | Merged | Nonchimera |
|---|---|---|---|---|---|---|---|
| 15E - 1 | AMB MP+ | 46968 | 44461 | 44392 | 44369 | 44242 | 44230 |
| 6B - 2 | OAW MP+ | 25343 | 23991 | 23882 | 23969 | 23811 | 23811 |
| 5B - 1 | OAW Control | 27398 | 26006 | 25982 | 25986 | 25976 | 25976 |
| 5aE - 1 | AMB MP+ | 45822 | 44152 | 44099 | 44109 | 43973 | 43973 |
| 10aA - 1 | OAW MP+ | 99390 | 94836 | 94653 | 94671 | 94353 | 94205 |
| 5aC - 2 | AMB Control | 77568 | 73174 | 73129 | 73130 | 73019 | 73019 |
| 15E - 2 | OAW Control | 23219 | 22069 | 22048 | 22045 | 22026 | 22026 |
| 15C - 1 | OAW MP+ | 82624 | 76841 | 76349 | 76438 | 73826 | 72033 |
| 10aB - 1 | AMB MP+ | 44719 | 41897 | 41807 | 41866 | 41495 | 41495 |
| 15E - 3 | AMB Control | 21016 | 19892 | 19813 | 19857 | 19742 | 19742 |
| 6D - 1 | AMB Control | 50868 | 48504 | 48363 | 48414 | 48090 | 48044 |
| 5aC - 1 | OAW Control | 33629 | 32048 | 31998 | 31994 | 31891 | 31859 |
| 11E - 2 | OAW MP+ | 93788 | 90475 | 90409 | 90402 | 90155 | 90146 |
| 5aA - 1 | OAW MP+ | 26111 | 24910 | 24846 | 24859 | 24708 | 24683 |
| 5D - 1 | AMB Control | 36048 | 34782 | 34765 | 34754 | 34710 | 34666 |
| 12B - 1 | OAW MP+ | 51802 | 50069 | 49983 | 50013 | 49874 | 49845 |
| 12E - 1 | AMB Control | 79689 | 76155 | 76112 | 76109 | 75821 | 75770 |
| 2E - 1 | AMB MP+ | 36469 | 34374 | 34349 | 34353 | 34289 | 34289 |
| 2B - 1 | AMB Control | 93851 | 90460 | 90421 | 90323 | 90201 | 90201 |
| 6B - 1 | AMB MP+ | 24767 | 23566 | 23506 | 23538 | 23364 | 23328 |
| 2aE - 1 | OAW Control | 56820 | 53651 | 53600 | 53596 | 53342 | 53342 |
| 10aC - 1 | OAW Control | 114866 | 110660 | 110347 | 110599 | 110242 | 110229 |
| 12E - 2 | AMB MP+ | 16088 | 14928 | 14822 | 14864 | 14767 | 14690 |
| 11E - 1 | OAW Control | 31587 | 30116 | 30045 | 30016 | 29943 | 29901 |
| 5C - 1 | AMB MP+ | 29051 | 27927 | 27900 | 27896 | 27851 | 27851 |
| 2D - 1 | OAW MP+ | 75311 | 73356 | 73247 | 73255 | 72971 | 72971 |
| 6E - 1 | OAW Control | 24571 | 23594 | 23583 | 23578 | 23571 | 23571 |
| 2aC - 1 | AMB MP+ | 21209 | 20343 | 20305 | 20329 | 20214 | 20185 |
| 10aD - 1 | AMB Control | 26776 | 25590 | 25551 | 25561 | 25491 | 25491 |
| 11A - 1 | AMB MP+ | 10752 | 10144 | 10095 | 10134 | 10026 | 10026 |
| 5C - 2 | OAW MP+ | 43598 | 42054 | 42021 | 42021 | 41807 | 41723 |
| 2A - 1 | OAW Control | 26542 | 25161 | 25122 | 25096 | 25002 | 25002 |
| 11B - 1 | AMB Control | 76572 | 73293 | 73263 | 73255 | 72994 | 72994 |
| 2aB - 1 | OAW MP+ | 38502 | 37175 | 37138 | 37084 | 36876 | 36876 |
| 12A - 1 | OAW Control | 31488 | 30216 | 30071 | 30042 | 29692 | 29692 |
| 2aD - 1 | AMB Control | 44886 | 43608 | 43583 | 43597 | 43543 | 43525 |
| NA | blank_control | 1918 | 1795 | 1792 | 1792 | 1792 | 1792 |
| NA | blank_control | 169 | 162 | 153 | 158 | 153 | 153 |
It is common at this point, especially in 16S/18S/ITS amplicon sequencing, to assign taxonomy to the sequence variants. The DADA2 package provides a native implementation of the naive Bayesian classifier method for this purpose. The assignTaxonomy function takes as input a set of sequences to be classified and a training set of reference sequences with known taxonomy, and outputs taxonomic assignments with at least minBoot bootstrap confidence.
The dada2 package GitHub maintains the most updated versions of the Silva databases., but I downloaded the databases from the associated Zenodo. The versions in this GitHub repository, used here, were last updated on 26 July 2022.
Let’s inspect the taxonomic assignments:
## Kingdom Phylum Class Order
## [1,] "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rickettsiales"
## [2,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Pseudomonadales"
## [3,] "Bacteria" "Cyanobacteria" "Cyanobacteriia" "Chloroplast"
## [4,] "Bacteria" "Myxococcota" "Myxococcia" "Myxococcales"
## [5,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Burkholderiales"
## [6,] "Bacteria" "Campylobacterota" "Campylobacteria" "Campylobacterales"
## Family Genus Species
## [1,] "Fokiniaceae" "MD3-55" NA
## [2,] "Pseudomonadaceae" "Pseudomonas" NA
## [3,] NA NA NA
## [4,] "Myxococcaceae" "P3OB-42" NA
## [5,] "Alcaligenaceae" "Alcaligenes" NA
## [6,] NA NA NA
Great! We can now save this and hand it off to Phyloseq for further analyses.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3277 taxa and 38 samples ]
## sample_data() Sample Data: [ 38 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 3277 taxa by 7 taxonomic ranks ]
Removal of mitochondira, chloroplasts, and non-bacteria taxa reduced the total number of taxa from 3493 to 3277.
Just 5 out of the 3277 ASVs were classified as contaminants.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3272 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 3272 taxa by 7 taxonomic ranks ]
Looks good! We have now cleaned up the sample data to remove contamination from non target organisms and those from the negative controls.
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3230 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 3230 taxa by 8 taxonomic ranks ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3110 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 3110 taxa by 8 taxonomic ranks ]
## [1] "samples with counts below z-score -2.5 :"
## character(0)
## [1] "zscores:"
## named numeric(0)
## [1] "OTUs passing frequency cutoff 1e-04 : 257"
## [1] "OTUs with counts in 0.02 of samples:"
##
## TRUE
## 257
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 257 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 257 taxa by 8 taxonomic ranks ]
## named numeric(0)
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 257 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 257 taxa by 8 taxonomic ranks ]
Notes from phyloseq author Visualize alpha-diversity - Should be done on raw, untrimmed dataset
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 8 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 5 sample variables ]
## tax_table() Taxonomy Table: [ 8 taxa by 8 taxonomic ranks ]
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 249 taxa and 36 samples ]
## sample_data() Sample Data: [ 36 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 249 taxa by 8 taxonomic ranks ]
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE
## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.002676 0.0008919 0.5372 0.6602
## Residuals 32 0.053128 0.0016602
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.002676 0.0008919 0.5372 999 0.666
## Residuals 32 0.053128 0.0016602
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## AMB Control AMB MP+ OAW Control OAW MP+
## AMB Control 0.67200 0.40300 0.674
## AMB MP+ 0.67533 0.29600 0.452
## OAW Control 0.39916 0.28025 0.718
## OAW MP+ 0.66196 0.45702 0.71966
##
## Call:
## adonis(formula = seq.acc ~ Treatment, data = samdf.acc, permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Treatment 3 1.1253 0.37509 0.92281 0.07962 0.746
## Residuals 32 13.0071 0.40647 0.92038
## Total 35 14.1323 1.00000
Taking phyloseq data and making some preliminary visualizations based on DADA2 tutorial:
## Run 0 stress 0.04934965
## Run 1 stress 0.0746166
## Run 2 stress 0.05637317
## Run 3 stress 0.05589854
## Run 4 stress 0.05199764
## Run 5 stress 0.07342775
## Run 6 stress 0.05192769
## Run 7 stress 0.05359359
## Run 8 stress 0.07439729
## Run 9 stress 0.05200148
## Run 10 stress 0.05352297
## Run 11 stress 0.05258621
## Run 12 stress 0.0739395
## Run 13 stress 0.04990088
## Run 14 stress 0.05273297
## Run 15 stress 0.07420514
## Run 16 stress 0.05318327
## Run 17 stress 0.07436481
## Run 18 stress 0.05389317
## Run 19 stress 0.05376737
## Run 20 stress 0.07443605
## Run 21 stress 0.0744205
## Run 22 stress 0.0523168
## Run 23 stress 0.05391028
## Run 24 stress 0.05120502
## Run 25 stress 0.05272165
## Run 26 stress 0.05487721
## Run 27 stress 0.05201142
## Run 28 stress 0.05261303
## Run 29 stress 0.05558435
## Run 30 stress 0.05181876
## Run 31 stress 0.05472496
## Run 32 stress 0.05175883
## Run 33 stress 0.05219736
## Run 34 stress 0.0742094
## Run 35 stress 0.07477284
## Run 36 stress 0.05243271
## Run 37 stress 0.07443597
## Run 38 stress 0.05534114
## Run 39 stress 0.05264528
## Run 40 stress 0.05476336
## Run 41 stress 0.05150825
## Run 42 stress 0.07425756
## Run 43 stress 0.05600592
## Run 44 stress 0.05493166
## Run 45 stress 0.07450414
## Run 46 stress 0.05180725
## Run 47 stress 0.05230675
## Run 48 stress 0.05396425
## Run 49 stress 0.05038315
## Run 50 stress 0.05254027
## Run 51 stress 0.07303486
## Run 52 stress 0.05685351
## Run 53 stress 0.05219459
## Run 54 stress 0.0517006
## Run 55 stress 0.05252119
## Run 56 stress 0.05259551
## Run 57 stress 0.05254038
## Run 58 stress 0.05280253
## Run 59 stress 0.07414479
## Run 60 stress 0.05386826
## Run 61 stress 0.05553755
## Run 62 stress 0.05533509
## Run 63 stress 0.07459378
## Run 64 stress 0.05169705
## Run 65 stress 0.05280704
## Run 66 stress 0.05050165
## Run 67 stress 0.05492977
## Run 68 stress 0.05366658
## Run 69 stress 0.05191697
## Run 70 stress 0.07335709
## Run 71 stress 0.07358913
## Run 72 stress 0.05201699
## Run 73 stress 0.07327029
## Run 74 stress 0.05192783
## Run 75 stress 0.07334119
## Run 76 stress 0.05178969
## Run 77 stress 0.052738
## Run 78 stress 0.05687152
## Run 79 stress 0.07441459
## Run 80 stress 0.0517904
## Run 81 stress 0.05174593
## Run 82 stress 0.07361126
## Run 83 stress 0.0744881
## Run 84 stress 0.0533005
## Run 85 stress 0.05217978
## Run 86 stress 0.05261949
## Run 87 stress 0.0741987
## Run 88 stress 0.05175881
## Run 89 stress 0.0739276
## Run 90 stress 0.05481606
## Run 91 stress 0.07430103
## Run 92 stress 0.0540038
## Run 93 stress 0.07475465
## Run 94 stress 0.0733594
## Run 95 stress 0.074414
## Run 96 stress 0.0549934
## Run 97 stress 0.05169708
## Run 98 stress 0.05329238
## Run 99 stress 0.05346095
## Run 100 stress 0.05691262
## Run 101 stress 0.05193697
## Run 102 stress 0.0517459
## Run 103 stress 0.05344187
## Run 104 stress 0.07380565
## Run 105 stress 0.05596173
## Run 106 stress 0.05278124
## Run 107 stress 0.05346922
## Run 108 stress 0.05174831
## Run 109 stress 0.04963303
## ... Procrustes: rmse 0.01670541 max resid 0.07686382
## Run 110 stress 0.07332924
## Run 111 stress 0.05503956
## Run 112 stress 0.07487197
## Run 113 stress 0.05516338
## Run 114 stress 0.05522886
## Run 115 stress 0.07473046
## Run 116 stress 0.05332271
## Run 117 stress 0.07454697
## Run 118 stress 0.05201137
## Run 119 stress 0.0740993
## Run 120 stress 0.07453052
## Run 121 stress 0.05199775
## Run 122 stress 0.05332381
## Run 123 stress 0.0537497
## Run 124 stress 0.05166873
## Run 125 stress 0.05043322
## Run 126 stress 0.07455454
## Run 127 stress 0.07455015
## Run 128 stress 0.07486829
## Run 129 stress 0.07421606
## Run 130 stress 0.05505515
## Run 131 stress 0.07447215
## Run 132 stress 0.05696038
## Run 133 stress 0.0745107
## Run 134 stress 0.05173746
## Run 135 stress 0.05242349
## Run 136 stress 0.05172763
## Run 137 stress 0.0738539
## Run 138 stress 0.05167663
## Run 139 stress 0.07452003
## Run 140 stress 0.0519278
## Run 141 stress 0.07545263
## Run 142 stress 0.07449724
## Run 143 stress 0.0739279
## Run 144 stress 0.05099502
## Run 145 stress 0.05180026
## Run 146 stress 0.07393434
## Run 147 stress 0.05451018
## Run 148 stress 0.05219447
## Run 149 stress 0.05158206
## Run 150 stress 0.05507125
## Run 151 stress 0.07441592
## Run 152 stress 0.07430414
## Run 153 stress 0.0549316
## Run 154 stress 0.07450411
## Run 155 stress 0.05256649
## Run 156 stress 0.05219474
## Run 157 stress 0.05096078
## Run 158 stress 0.05352463
## Run 159 stress 0.05483983
## Run 160 stress 0.05230158
## Run 161 stress 0.07449742
## Run 162 stress 0.05270281
## Run 163 stress 0.05192775
## Run 164 stress 0.05144565
## Run 165 stress 0.07351521
## Run 166 stress 0.05222437
## Run 167 stress 0.05492416
## Run 168 stress 0.0550851
## Run 169 stress 0.05641403
## Run 170 stress 0.05238796
## Run 171 stress 0.05192782
## Run 172 stress 0.07458608
## Run 173 stress 0.05443599
## Run 174 stress 0.05230684
## Run 175 stress 0.05553766
## Run 176 stress 0.07357158
## Run 177 stress 0.0519976
## Run 178 stress 0.05509418
## Run 179 stress 0.07430424
## Run 180 stress 0.05172717
## Run 181 stress 0.05273072
## Run 182 stress 0.05112478
## Run 183 stress 0.05420418
## Run 184 stress 0.05372851
## Run 185 stress 0.05041913
## Run 186 stress 0.05227609
## Run 187 stress 0.05493204
## Run 188 stress 0.05256644
## Run 189 stress 0.0493495
## ... New best solution
## ... Procrustes: rmse 8.487414e-05 max resid 0.0002225465
## ... Similar to previous best
## Run 190 stress 0.05258621
## Run 191 stress 0.05228849
## Run 192 stress 0.07475731
## *** Solution reached
Figure SXX. Relative abundance of major ITS2 types by coral fragment (x axis) and treatment (facets). Light green represent Symbiodinium spp. (A3) and dark green represents Breviolum spp. (B2).
Figure SXX. Relative abundance of major ITS2 types grouped by treatment (n = 4-5 corals per treatment). Light green represent Symbiodinium spp. (A3) and dark green represents Breviolum spp. (B2).
All code was written by Colleen B. Bove, feel free to contact with questions.
Session information from the last run date on 29 July 2022:
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] awtools_0.2.1 microbiome_1.8.0
## [3] ggpubr_0.4.0 MCMC.OTU_1.0.10
## [5] MCMCglmm_2.33 ape_5.6-1
## [7] coda_0.19-4 Matrix_1.3-4
## [9] janitor_2.1.0 vegan_2.5-7
## [11] lattice_0.20-45 permute_0.9-7
## [13] plotly_4.10.0 RColorBrewer_1.1-3
## [15] decontam_1.6.0 phyloseq_1.30.0
## [17] kableExtra_1.3.4 ShortRead_1.44.3
## [19] GenomicAlignments_1.22.1 SummarizedExperiment_1.16.1
## [21] DelayedArray_0.12.3 matrixStats_0.62.0
## [23] Biobase_2.46.0 Rsamtools_2.2.3
## [25] GenomicRanges_1.38.0 GenomeInfoDb_1.22.1
## [27] Biostrings_2.54.0 XVector_0.26.0
## [29] IRanges_2.20.2 S4Vectors_0.24.4
## [31] BiocParallel_1.20.1 BiocGenerics_0.32.0
## [33] forcats_0.5.1 stringr_1.4.0
## [35] dplyr_1.0.8 purrr_0.3.4
## [37] readr_2.1.2 tidyr_1.2.0
## [39] tibble_3.1.7 ggplot2_3.3.6
## [41] tidyverse_1.3.1 dada2_1.20.0
## [43] Rcpp_1.0.9 knitr_1.33
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.4.1 systemfonts_1.0.2
## [4] plyr_1.8.7 igraph_1.2.11 lazyeval_0.2.2
## [7] splines_3.6.3 crosstalk_1.2.0 usethis_2.0.1
## [10] digest_0.6.29 foreach_1.5.2 htmltools_0.5.2
## [13] fansi_1.0.3 memoise_2.0.1 magrittr_2.0.3
## [16] cluster_2.1.2 remotes_2.4.2 tzdb_0.2.0
## [19] modelr_0.1.8 RcppParallel_5.1.5 svglite_2.1.0
## [22] prettyunits_1.1.1 jpeg_0.1-9 colorspace_2.0-3
## [25] rvest_1.0.1 textshaping_0.3.6 haven_2.4.3
## [28] xfun_0.29 callr_3.7.0 crayon_1.5.1
## [31] RCurl_1.98-1.7 jsonlite_1.7.2 survival_3.2-13
## [34] iterators_1.0.14 glue_1.6.2 gtable_0.3.0
## [37] zlibbioc_1.32.0 webshot_0.5.2 pkgbuild_1.3.1
## [40] car_3.1-0 Rhdf5lib_1.8.0 abind_1.4-5
## [43] scales_1.2.0 DBI_1.1.3 rstatix_0.7.0
## [46] viridisLite_0.4.0 htmlwidgets_1.5.4 httr_1.4.3
## [49] ellipsis_0.3.2 farver_2.1.1 pkgconfig_2.0.3
## [52] sass_0.4.0 dbplyr_2.1.1 deldir_1.0-6
## [55] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1
## [58] rlang_1.0.4 reshape2_1.4.4 munsell_0.5.0
## [61] cellranger_1.1.0 tools_3.6.3 cachem_1.0.6
## [64] cli_3.3.0 generics_0.1.3 ade4_1.7-18
## [67] devtools_2.4.2 broom_1.0.0 evaluate_0.15
## [70] biomformat_1.14.0 fastmap_1.1.0 ragg_1.1.3
## [73] yaml_2.3.5 processx_3.5.2 fs_1.5.2
## [76] nlme_3.1-155 xml2_1.3.3 brio_1.1.3
## [79] compiler_3.6.3 rstudioapi_0.13 curl_4.3.2
## [82] png_0.1-7 testthat_3.1.2 ggsignif_0.6.3
## [85] reprex_2.0.1 bslib_0.4.0 stringi_1.7.8
## [88] ps_1.6.0 highr_0.9 desc_1.4.1
## [91] cubature_2.0.4.2 tensorA_0.36.2 multtest_2.42.0
## [94] vctrs_0.4.1 pillar_1.8.0 lifecycle_1.0.1
## [97] jquerylib_0.1.4 cowplot_1.1.1 data.table_1.14.2
## [100] bitops_1.0-7 corpcor_1.6.10 R6_2.5.1
## [103] latticeExtra_0.6-30 hwriter_1.3.2.1 sessioninfo_1.2.2
## [106] codetools_0.2-18 pkgload_1.2.1 MASS_7.3-55
## [109] assertthat_0.2.1 rhdf5_2.30.1 rprojroot_2.0.3
## [112] withr_2.5.0 GenomeInfoDbData_1.2.2 mgcv_1.8-38
## [115] hms_1.1.1 grid_3.6.3 rmarkdown_2.13
## [118] snakecase_0.11.0 carData_3.0-5 Rtsne_0.15
## [121] lubridate_1.8.0 interp_1.0-33